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In  the  current  work,  a  computational  model  of  a  microfluidic  fuel  cell  with  flow-through  porous  elec¬ 
trodes  is  developed  and  validated  with  experimental  data  based  on  vanadium  redox  electrolyte  as  fuel 
and  oxidant.  The  model  is  the  first  of  its  kind  for  this  innovative  fuel  cell  design.  The  coupled  prob¬ 
lem  of  fluid  flow,  mass  transport  and  electrochemical  kinetics  is  solved  from  first  principles  using  a 
commercial  multiphysics  code.  The  performance  characteristics  of  the  fuel  cell  based  on  polarization 
curves,  single  pass  efficiency,  fuel  utilization  and  power  density  are  predicted  and  theoretical  maxima 
are  established.  Fuel  and  oxidant  flow  rate  and  its  effect  on  cell  performance  is  considered  and  an  optimal 
operating  point  with  respect  to  both  efficiency  and  power  output  is  identified  for  a  given  flow  rate.  The 
results  help  elucidate  the  interplay  of  kinetics  and  mass  transport  effects  in  influencing  porous  electrode 
polarization  characteristics.  The  performance  and  electrode  polarization  at  the  mass  transfer  limit  are 
also  detailed.  The  results  form  a  basis  for  determining  parameter  variations  and  design  modifications  to 
improve  performance  and  fuel  utilization.  The  validated  model  is  expected  to  become  a  useful  design  tool 
for  development  and  optimization  of  fuel  cells  and  electrochemical  sensors  incorporating  microfluidic 
flow-through  porous  electrodes. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Microscale  power  sources  are  currently  under  rapid  develop¬ 
ment  due  to  the  high  demand  for  use  in  small,  portable  applications. 
These  applications  include  personal  electronic  devices  such  as 
mobile  phones,  laptops  as  well  as  wireless  sensors,  global  posi¬ 
tioning  systems,  medical  diagnostics  and  other  specialized  devices. 
Next  generation  power  sources  serving  these  applications  are 
required  to  be  light  weight,  have  high  power  densities,  high  durabil¬ 
ity  and  long  periods  of  operation  between  recharge.  Miniature  fuel 
cells  [1-3]  have  shown  great  promise  as  power  sources  for  small- 
scale  applications  due  to  their  inherently  high  power  and  energy 
density.  However,  scaling  down  conventional  fuel  cell  technolo¬ 
gies  to  produce  miniature  power  sources  has  several  challenges. 
Mechanical  constraints  arise  due  to  limitations  in  graphite  plate 
size  that  can  be  machined,  and  miniaturization  decreases  the  struc¬ 
tural  strength  of  the  membrane  electrode  assembly.  Furthermore, 
hydrogen  based  fuel  cells  require  hydrogen  storage  units  that  are 
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generally  too  bulky  for  portable  applications.  Direct  liquid  fuel 
cells  on  the  other  hand  suffer  low  cell  voltage  due  to  slow  elec¬ 
trochemical  kinetics  and  fuel  crossover.  There  are  also  noteworthy 
challenges  associated  with  membrane  hydration  and  degradation 
common  to  all  polymer  electrolyte  membrane  fuel  cell  technolo¬ 
gies. 

A  recent,  novel  fuel  cell  architecture  based  on  microfluidic  lab- 
on-chip  devices  has  shown  great  promise  in  addressing  many  of  the 
above  problems.  These  microfluidic  or  laminar  flow  based  fuel  cells 
incorporate  all  the  fundamental  components  of  a  conventional  fuel 
cell  to  a  microfluidic  channel  and  its  walls  [4].  Such  fuel  cells  oper¬ 
ate  without  a  membrane  and  in  the  most  common  configuration 
use  the  laminar  nature  of  microscale  flows  to  maintain  separation 
between  fuel  and  oxidant  streams,  and  at  the  same  time  facilitate 
ionic  charge  transfer.  Crossover  effects  are  eliminated  by  strategic 
cell  design  provided  that  interdiffusion  is  restricted  to  a  small  inter¬ 
facial  width  at  the  center  channel,  from  which  the  electrodes  are 
positioned  sufficiently  far  away.  In  this  way,  microfluidic  fuel  cells 
prove  to  have  a  number  of  unique  advantages  [4,5]:  the  need  for  an 
ion-exchange  membrane  is  eliminated  and  along  with  it  the  mem¬ 
brane  related  issues;  sealing,  manifolding  and  fuel  delivery  system 
requirements  are  reduced;  manufacturing  of  these  fuel  cells  using 
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Fig.  1.  Schematic  of  the  microfluidic  fuel  cell  with  flow-through  porous  electrodes,  showing  the  3D  and  2D  modeling  domains  in  parts  (a)  and  (b),  respectively. 


established  microfabrication  methods  is  inexpensive  [6];  and  they 
can  operate  at  room  temperature  and  do  not  require  any  auxiliary 
systems  for  water  management  and  cooling. 

Notwithstanding  these  advantages,  as  with  any  emerging  tech¬ 
nology,  many  challenges  and  opportunities  for  improvement  exist. 
The  lack  of  fuel  distribution  is  an  important  consideration  for 
microfluidic  fuel  cells  as  well  as  other  miniaturized  fuel  cells,  and 
other  issues  such  as  sealing  arise  due  to  use  of  liquid  fuels  and  elec¬ 
trolytes.  However,  microfluidic  fuel  cells  with  flow-through  porous 
electrodes  can  conveniently  be  recharged  in  situ  by  running  the 
cell  in  electrolytic  mode  [7]  in  a  closed-loop  system,  thus  eliminat¬ 
ing  the  need  for  frequent  refueling.  Moreover,  although  recharging 
is  expected  to  be  slower  than  manually  replacing  a  fuel  cartridge, 
it  can  potentially  be  faster  than  for  conventional  batteries  due  to 
the  convective  mass  transport  of  the  flow-through  cell  architec¬ 
ture.  On  the  system  level,  microfluidic  fuel  cells  require  a  pump  to 
operate.  While  vanadium  redox  based  electrochemical  cells  have 
relatively  high  durability  [8],  the  pump  may  have  a  negative  effect 
on  overall  fuel  cell  system  lifetime  and  may  need  to  be  replaced 
more  frequently  than  the  electrolyte.  The  cost  and  parasitic  power 
requirements  of  the  pump  are  also  important  to  consider  for  sys¬ 
tem  integration  of  microfluidic  fuel  cells.  Other  challenges  include 
high  contact  resistances  which  result  from  establishing  electrical 
contact  between  porous  electrodes  and  the  external  circuit. 

A  variety  of  microfluidic  fuel  cells  have  been  demonstrated 
using  different  fuel  and  oxidant  combinations,  electrode  materials 
and  channel  configurations  [4].  One  of  the  first  such  devices  was 
reported  by  Ferrigno  et  al.  [9],  utilizing  vanadium  species  for  both 
half-cells  in  a  Y-shaped  microchannel.  Using  a  similar  Y-shaped 
channel  design,  Choban  et  al.  [10]  demonstrated  a  microfluidic  fuel 
cell  running  on  formic  acid  as  fuel.  Cohen  et  al.  [1 1  ]  introduced  the 
planar  or  F-shaped  microchannel  design  fuel  cell  using  formic  acid 
as  fuel.  This  design  had  increased  contact  area  between  the  elec¬ 
trolyte  and  electrodes  as  compared  to  previous  designs.  Devices 
using  oxygen  as  the  oxidant  suffered  from  low  power  densities 
due  to  the  poor  solubility  of  oxygen  in  liquid  electrolytes.  This 
problem  was  addressed  by  Jayashree  et  al.  [12],  who  introduced 
the  novel,  air  breathing  electrode  architecture.  All  the  above  men¬ 
tioned  devices  were  reported  to  be  limited  by  mass  transport  to 
the  active  sites,  and  as  a  result  had  low  fuel  utilization  and  power 
density  values.  Kjeang  et  al.  [7,13,14]  recently  reported  incremental 
design  improvements  to  mitigate  this  limitation,  including  graphite 
rod  arrays  and  porous  carbon  paper  as  electrodes,  culminating  in  a 
microfluidic  fuel  cell  with  flow-through  porous  carbon  electrodes 
using  vanadium  species  in  both  half-cells.  This  device  has  demon¬ 
strated  some  of  the  highest  power  density  levels  (131  mWcrn-2) 


and  fuel  utilization  values  (100%)  of  ah  the  devices  reported  till 
date  [3-5]  along  with  in  situ  fuel  regeneration  capabilities.  More 
recently,  a  few  other  architectures  have  been  reported  by  Salloum 
et  al.  [  1 5-1 7]  using  flow-through  electrodes  and  vanadium  species. 

Modeling  of  microfluidic  fuel  cells  is  important  to  cut  down  the 
time  involved  in  prototyping,  building  and  characterizing  actual 
devices,  and  to  understand  the  effect  of  various  design  parame¬ 
ters  on  the  performance.  Modeling  can  also  elucidate  the  complex 
mass  transfer  and  electrochemical  coupling  effects  which  exist  in 
these  devices.  The  first  Computational  Fluid  Dynamics  (CFD)  based 
model  for  microfluidic  fuel  cells  was  developed  by  Bazylak  et  al.  [  1 8  ] 
who  considered  a  T-shaped  microchannel  and  suggested  methods 
to  improve  fuel  utilization  by  using  tapered  electrodes.  Chang  et  al. 
[19]  presented  a  theoretical  model  which  included  a  Butler-Volmer 
model  for  the  electrochemical  kinetics  in  a  Y-shaped  microchannel. 
A  similar  model  for  the  planar  or  F-shaped  channel  was  reported  by 
Chen  et  al.  [20]  and  validated  based  on  the  experiments  performed 
by  Cohen  et  al.  [  1 1  ].  No  theoretical/numerical  models  have  yet  been 
developed  for  the  flow-through  electrode  architecture. 

The  overall  objective  of  the  present  work  is  to  develop  and  val¬ 
idate  a  detailed  computational  model  for  a  microfluidic  fuel  cell 
with  flow-through  porous  electrodes.  The  model  incorporates  ah 
the  fundamental  phenomena  in  the  fuel  cell  including  fluid  flow  in 
microchannels  and  porous  media,  mass  transport  and  electrochem¬ 
ical  kinetics,  and  resolves  the  characteristics  of  both  individual 
half-cells  (anode  and  cathode)  as  well  as  the  combined  fuel  cell. 
Such  a  model  is  expected  to  be  a  vital  tool  to  understand  the 
unique  physicochemical  characteristics  of  this  configuration  and 
to  support  the  development  and  optimization  of  cell  architectures. 
Results  are  expected  to  shed  light  on  electrode  polarization  char¬ 
acteristics  and  guide  the  design  of  higher  performance  electrode 
architectures.  In  addition,  this  work  is  intended  to  provide  a  gen¬ 
eral  numerical  modeling  framework  for  electrochemical  devices 
with  microfluidic  flow-through  architectures,  such  as  microfluidic 
bio-sensors. 

2.  Model  formulation  and  theory 

The  microfluidic  fuel  cell  architecture  considered  in  this  work  is 
based  on  vanadium  redox  reactions  at  both  half-cells  and  incor¬ 
porates  porous  flow-through  carbon  paper  electrodes  [7].  A  full 
numerical  description  of  all  fundamental  transport  and  reaction 
phenomena  occurring  in  the  fuel  cell  is  highly  complex  with  a  large 
degree  of  non-linearity,  and  to  obtain  a  solution  it  is  necessary  to 
simplify  the  problem  to  a  manageable  level.  In  order  to  present  a 
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straightforward  yet  accurate  model,  the  main  assumptions  made 
are  listed  below: 

1.  Isothermal  and  steady  state  conditions  are  considered. 

2.  Incompressible  fluid  flow  is  assumed  and  effects  of  gravity  are 
neglected. 

3.  Physical  properties  of  the  electrode  and  electrolyte  are  consid¬ 
ered  homogeneous. 

4.  Dilute  solution  approximation  is  assumed  to  hold. 

5.  Ionic  migration  is  neglected  owing  to  the  high  concentrations  of 
supporting  electrolyte. 

Assumptions  4  and  5  are  justified  owing  to  the  aqueous,  acidic 
nature  of  the  electrolyte  used  in  [7]  and  modeled  here.  Apart  from 
these  general  assumptions,  others  regarding  the  electrochemical 
kinetics  and  mass  transport  aspects  of  the  model  will  be  detailed 
in  following  sections  as  they  arise. 

2.1.  Geometry 

A  schematic  of  the  fuel  cell  indicating  the  modeling  domain  is 
shown  in  Fig.  1.  The  electrolyte  streams  enter  the  fuel  cell  via  the 
two  inlets  on  either  side,  and  after  passing  through  the  inlet  reser¬ 
voirs,  flow  through  the  porous  electrodes.  The  two  streams  exit 
the  electrodes  into  the  common  center  channel  where  they  make 
a  90°  turn  and  flow  in  a  stratified,  co-laminar  format  towards  the 
outlet.  The  channel  thickness  is  constant  at  300  p,m.  Due  to  plug 
flow  conditions  in  the  porous  electrodes,  it  is  deemed  sufficient  to 
model  a  two  dimensional  (2D)  slice  for  the  mass  transport  and  elec¬ 
trochemistry.  However,  variations  in  flow  velocity  at  the  entrance 
to  the  electrodes  along  their  length  (y  direction)  are  expected.  For 
effectively  capturing  the  flow  distribution,  the  fluid  flow  is  solved 
using  a  full  three  dimensional  (3D)  model  of  the  cell.  The  resulting 
3D  calculations  are  used  to  define  the  fluid  flow  boundary  condi¬ 
tions  at  the  inlet  of  the  electrodes  in  the  2D  model.  The  domain  used 
for  3D  modeling  is  shown  in  Fig.  1  (a)  and  the  same  for  the  2D  model 
in  Fig.  1(b).  Note  that  the  cell  is  bisected  along  the  channel  thick¬ 
ness  for  the  3D  model,  to  exploit  the  symmetry  about  this  plane. 
The  position  of  the  current  collectors  at  the  end  of  each  electrode 
is  shown  in  Fig.  1.  The  mid-plane  bisecting  the  cell  width  wise  is 
used  for  half-cell  modeling  purposes. 

2.2.  Fluid  flow  model 

A  complete  3D,  laminar  CFD  model  is  used  to  capture  the  flow 
distribution  in  the  fuel  cell.  The  carbon  paper  electrodes  are  treated 
as  porous  bodies  with  homogeneous  permeability  in  the  in-plane 
direction.  The  pressure  drop  across  the  electrodes  is  calculated 
using  the  well-known  Darcy’s  law  [21  ].  The  continuity  and  momen¬ 
tum  equations,  written  for  regions  inside  the  electrodes  are: 

V-i)  =  0  (1) 

Vp  =  -£t)  (2) 

where  v  is  the  superficial  fluid  velocity  vector,  p  is  the  pressure  in 
the  liquid,  fi  is  the  dynamic  viscosity  and  I<  is  the  permeability  of  the 
carbon  paper  in  the  in-plane  direction  [22,23].  The  incompressible 
equations  for  continuity  and  momentum  conservation  [21]  used 
for  the  flow  in  the  regions  outside  the  porous  electrodes  are  given 


by: 

V.v  =  0 

(3) 

dv  ^ .  1  Li 

T-  +  (v.Vv)  =  --S/p  +  —Vzv 

3 1  p  p 

(4) 

2.3.  Electrochemistry  and  mass  transport  model 


The  governing  equations  for  species  conservation  and  charge 
conservation  in  the  carbon  paper  electrode  and  electrolyte  region 

are  written  as  [24,25]: 

vVCj-Dfv2Cj  =  Sj 

(5) 

Vis  =  -Vi,  =  -afw2<t>s  =  kfv 2<t>, 

(6) 

The  variables  in  Eqs.  (5)  and  (6)  are  explained  as  follows:  Cj  denotes 
the  bulk  concentration  of  the  vanadium  species  j,  where  j  e  {II,  III ,  IV 
and  V}.  The  symbols  {II,  III,  IV  and  V}  refer  to  the  vanadium  species 
{ V2+ ,  V3+ ,  V02+  and  VO +}  respectively,  in  all  the  following  sec¬ 
tions.  Djff  are  the  effective  diffusion  coefficients  of  species  j;  Sj  are 
the  species  source  terms  for  the  species  j;  <ps,  0/  are  the  potentials 
in  the  electrode  and  electrolyte,  respectively;  h  are  the  electrode 
and  electrolyte  current  densities;  and  a and  kf^  are  the  effective 
conductivities  of  the  electrode  and  electrolyte,  respectively. 

Mixed  evidence  [26-28]  exists  regarding  redox  kinetics  of  vana¬ 
dium  reactions  occurring  on  carbon  electrodes.  These  have  been 
characterized  to  have  a  complex  mechanism,  with  combined  chem¬ 
ical  and  electron  transfer  reactions  [26-28],  of  which  the  single 
electron  transfer  is  the  rate  determining  step  [26].  Moreover,  the 
exact  nature  of  the  solvated  vanadium  ions  has  been  reported  to 
be  dependent  on  pH  and  vanadium  ion  concentration  [29].  How¬ 
ever,  a  simplified  treatment  of  the  redox  kinetics  involving  single 
electron  transfer  reactions  has  been  demonstrated  [25]  to  be  suffi¬ 
cient  in  capturing  the  qualitative  cell  behaviour.  The  overall  redox 
reactions  considered  here  are: 

VO+ +  2H+ +  e~  ^±V02+ +  H20,  11°  =  0.750  V  vs  SCE  (7) 

V3+  +  e~  -  V2+ ,  D°  =  -0.496  V  vs  SCE  (8) 

where  D°  and  U°  are  the  standard  equilibrium  potentials  versus 
the  Saturated  Calomel  Electrode  (SCE)  [30].  Other  than  the  primary 
reactions  given  above,  a  secondary  reaction  involving  the  V3+  and 
V02+  species  is  known  to  occur  [31  ]: 

V02+ +  2H+ +  e~  ^V3+ +H20,  U°  =  0.07  V  vs  SCE  (9) 

During  normal  operating  cell  conditions,  this  reaction  has  a  neg¬ 
ligible  rate  [32  ]  due  to  slow  kinetics  and  low  reactant  concentration. 
Hydrogen  and  oxygen  evolution  may  also  occur  as  side  reactions  in 
the  acidic,  aqueous  electrolyte  present  in  the  cell.  However,  previ¬ 
ous  studies  show  that  the  kinetics  of  these  side  reactions  is  around 
five  orders  of  magnitude  lower  than  the  main  vanadium  redox 
reactions  (Eqs.  (7)  and  (8))  on  carbon  electrodes  [26,33,34].  Further¬ 
more,  during  typical  fuel  cell  operation  the  local  overpotentials  are 
not  high  enough  for  these  reactions  to  occur.  Therefore  from  ther¬ 
modynamic  and  kinetics  considerations  at  leading  order  the  effect 
of  these  reactions  is  negligible  and  consequently  not  included  in 
the  model. 

Experimental  investigations  by  Fabjan  et  al.  [26]  into  the  vana¬ 
dium  redox  kinetics  have  shown  that  a  single  electron  transfer  is 
the  rate  determining  step,  making  it  amenable  to  a  Butler-Volmer 
[35,36]  model.  Other  works  [27,28]  suggest  an  alternative  mech¬ 
anism  for  the  V02+/V0j  couple.  Several  earlier  modeling  efforts 
of  vanadium  redox  batteries  [25,37]  have  shown  that  performance 
of  a  cell  based  on  these  reactions  is  adequately  captured  by  the 
Butler-Volmer  equation.  Therefore,  given  the  uncertainties  that 
exist  with  regards  to  the  kinetics  of  the  vanadium  redox  reactions, 
the  Butler-Volmer  equation  [35,36]  is  deemed  a  necessary  and 
justified  simplification  for  the  present  model.  The  charge  transfer 
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current  densities  corresponding  to  the  redox  reactions  are  thus 
given  by: 


:  _  :0 
H  —  h 


:  _  :0 
12  —  ^2 


-«2,c^2 

RT 


(10) 


(11) 


Note  that  subscript  1  refers  to  the  V02+/V0j  side,  which  is  the 
cathode  during  fuel  cell  operation,  and  subscript  2  refers  to  the 
V2+IV3+  side,  or  the  anode.  The  superscript  s  refers  to  the  surface 
concentration  of  the  vanadium  species,  a  is  the  charge  transfer  coef¬ 
ficient,  with  subscripts  a  and  c  referring  to  oxidation  (anodic)  and 
reduction  (cathodic)  reactions,  respectively.  F is  Faraday’s  constant, 
R  is  the  universal  gas  constant,  T  is  the  operating  temperature  and 
rj  is  the  overpotential.  The  exchange  current  density  i°  represents 
the  rate  of  charge  transfer  when  the  reaction  is  at  equilibrium.  For 
the  cathode  and  anode  side  it  is  written  as: 


i°=Fk,(clvr-‘{cvr* 

%  =  Fk2(cnf*-'{cmf*-' 


(12) 

(13) 


where  k\  and  k2  are  the  standard  reaction  rate  constants.  The  over¬ 
potentials  at  the  cathode  and  anode  side  are  written  as: 

>h  =  4>s  -4>i-u,  (14) 

m  =  -4>i-u2  (15) 

The  Nernst  equation  is  used  to  estimate  the  half-cell  potentials: 

u'-u°+ ?'"(£)  (,6) 


RT  , 


U2  =  U“+—ln  f- 


(?) 


(17) 


ii  and  i2  are  used  to  define  the  source  terms  for  the  charge  con¬ 
servation  equations  (Eq.  (6))  in  the  electrode  and  electrolyte  which 
are  given  by: 


Vis  =  -V.q  =  S 1  =  -aii 
Vis  =  —  V.q  =S2  =  -ai2 


(18) 

(19) 


where  a  is  the  internal  active  surface  area  of  the  electrode,  the 
value  of  which  is  estimated  based  on  fibre  diameter  and  poros¬ 
ity  and  is  listed  in  Table  3.  The  vanadium  species  are  transported 
to  and  from  the  electrode  surface  through  the  combined  effect  of 


coefficients  of  the  species.  These  are  estimated  using  an  empirical 
relationship  for  mass  transfer  in  carbon  fibre  electrodes  [38]: 


Sh  =  7 Re0  4,  (56  <  Pe  <  280) 


(21) 


The  Sherwood  number  and  Reynolds  number  are  calculated 
with  respect  to  the  carbon  fibre  diameter  (d)  and  the  apparent 
(superficial)  velocity  in  the  porous  medium.  The  Peclet  number 
range  used  for  this  study  is  of  the  same  order  of  magnitude  as 
the  above  range,  making  this  a  suitable  correlation.  In  the  fluidic 
domain  inside  the  electrodes,  these  three  dimensionless  quantities 
are  written  as: 


Re  = 


P\v\d 

P 


Shi  =  ^ 


Pej  = 


\v\d 

D, 


(22) 

(23) 

(24) 


Expanding  Eq.  (21)  for  each  vanadium  species  using  Eqs. 
(22)-(24)  we  can  write, 


I  _7Di  (p\v\d\ 
m]  ~  d  l  ft  ) 


0.4 


(25) 


This  correlation  is  valid  for  mass  transfer  for  a  single  fibre  and  pro¬ 
vides  an  upper  limit  for  a  bundle  of  fibres.  The  agreement  with 
experimental  data  is  reported  to  be  within  30%  error  [38]. 

Since  the  redox  reactions  being  modeled  are  surface  based,  there 
is  a  direct  coupling  between  the  mass  transport  and  the  electro¬ 
chemistry.  At  steady  state,  the  rate  at  which  a  chemical  species  is 
transported  to  or  from  the  surface  of  the  carbon  fibres  equals  the 
rate  of  the  reaction.  Hence,  the  relationship  can  be  expressed  math¬ 
ematically  between  the  mass  transport  and  electrochemical  source 
terms.  From  Eqs.  (18)-(20): 

S,v  =  -Sv  =  y  (26) 

S„  =  -S„,  =  Sf  (27) 

Eqs.  (26)  and  (27)  also  express  the  relationship  between  surface  and 
bulk  concentration  for  the  reacting  species,  and  the  correspond¬ 
ing  overpotential  on  the  anode  and  cathode  side.  Rearranging  this 
equation,  we  can  write  the  surface  concentration  for  each  species 
as  a  function  of  the  bulk  concentrations  and  overpotentials  [25]: 


rs  - 
Liv  ~ 


CV  - 


CII  - 


rs  - 
liii  ~ 


(ki/femj)(c/vf,'c(cv),““,'“  exp((-a1,cF^1)/RT)  +  {l+(k1/kmj)(c,vf'.c(cvr“’.“  exp ((-ahcFm)/ RT)}clv 
1  +(ki/kmjXc/vr“l  c(cv)“’-“  exp((a1,aF^1)/RT)  +  (k1/kmj)(c/v)“'^(cv )-“’•“  exp((-a1,cF^1)/RT)  ’ 

(ki/kmj)(c/v)1-“1’c(cv)“1’“  exp ((-ahaFm)/RT)  +  {1  +  (k,/km])(clvra^(cvr’a exp ((a^Fij, )/RT)}cv 


j  e  IV 


1  +  (ki/M(c/vr“v(cv)“'-“  exp((alaF^1  )/RT)  +  (k1/kmj)(c,vf  ^(cv  )““’■“ 
(k2/kmjXc„r2-c(cu,y-a^  exp ((-a2,cF^2)/RT)  +  {1  +  (fc2/kmj)(c,,)“^(c„,r^. 


exp((— ai  ,cFih)/RT)  ’ 
exp((-a2xFri2)/RT)\Cii 


1  +  (fe /M(c«r“2’c(c///)“2’a  exp((a2,aF^2)/RT)  +  (k2/kmjXc„r2-c(c,„ra2’a  exp((-a2,cF/?2)/RT) 

(l<2/kmjXc,iy-a2^clllr^  exp ((a2,aF^2)/RT)  +  {1  +  (h/k^Xc,,)-012^,,,)01^  exp ((a2,„F^2)/RT)}cffl 
1  +(k2/kmjXc„ra2’c(c,„r2’a  exp((o!2iaF?)2)/RT)  +  (k2/kmj)(c//)“2'c(cjj/r“2-a  exp((-a2,cF/?2)/RT)  ’ 


jeV 


jell 


Jem 


(28) 


(29) 


(30) 


(31) 


convection  and  diffusion  inside  the  pores  of  the  electrodes.  This 
rate  of  transport  is  proportional  to  the  difference  in  concentration 
between  bulk  and  the  surface  and  defines  the  species  source  terms 
in  Eq.  (5)  as: 

Sj  =  akmj(Cj  —  Cj),  je  {II,  III,  IV  and  V}  (20) 

where  /<mj  are  the  mass  transfer  coefficients  for  each  vana¬ 
dium  species,  which  depend  on  the  flow  velocity  and  diffusion 


Several  works  [23,39-43]  have  determined  the  mass  transport 
related  effective  properties  of  porous  media.  For  TORAY  carbon 
paper,  studies  by  Zamel  et  al.  [39]  show  that  the  effective  in-plane 
diffusivity  predicted  by  Bruggeman’s  correction  [43]  is  within  15% 
of  different  experimental  and  numerical  studies  [23,39-41  ].  Inside 
the  carbon  paper  the  effective  diffusion  coefficients  are  calculated 
based  on  this  correction  [43]: 

De/f  =  D.-fi1  '5  (32) 
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where  s  is  the  porosity  of  the  carbon  paper.  The  electrolyte  con¬ 
ductivity  in  the  microchannels  and  electrode  pores  are  calculated 


as: 

k‘  =  ^EzfDJcJ 

j 

(33) 

r-wZW* 

(34) 

j 


Here,  je{II,  III,  IV,  V  and  H+}  to  account  for  the  electrolyte  con¬ 
ductivity  of  the  acid.  kt  is  the  electrolyte  conductivity  outside  the 
electrode  regions  and  Zj  is  the  charge  number  for  each  of  the  ionic 
species.  The  proton  concentration  in  the  solution  is  assumed  to  be 
constant  at  8  M  since  changes  in  proton  concentration  are  small 
during  fuel  cell  operation. 

2.4.  Boundary  conditions  and  parameters 


2.4. 1.  Fluid  flow  boundary  conditions 

In  the  3D  fluid  flow  model,  the  following  boundary  conditions 
are  used.  Since  the  geometry  exhibits  symmetry  about  the  half  z- 
plane,  we  can  reduce  the  model  domain  by  half  by  implementing  a 
symmetry  boundary  condition  at  this  plane  (Fig.  1(a)): 

^  =  0,  ^  =  0  at  z  =  0  for  all  x  and  y  (35) 

dz  dz 

At  the  inlet,  the  inward  normal  velocity  is  written  as  a  function 
of  volume  flow  rate  and  inlet  channel  area: 


VP 

in 

ain  ’ 

x  =  0, 

yi  <  y  <  l, 

n  h 

°<z<- 

V? 

in 

ain  ’ 

x  =  W, 

y^<y<L, 

0  <Z  <  - 

-  -  2 

where  V°  is  the  inlet  volume  flow  rate  for  the  cathode  and  anode 
side,  ain  is  the  area  of  the  inlet  channel  and  h  is  the  channel  thick¬ 
ness.  The  outlet  is  specified  as  a  constant  pressure  boundary: 

p  =  ojx2<x<x3,  y  =  0,  0  <  z  <  ^  |  (37) 

For  the  2D  model,  the  flow  field  results  from  the  3D  model  are 
averaged  across  the  channel  thickness  and  mapped  as  a  bound¬ 
ary  condition  at  the  inlet  to  the  electrodes.  Expressing  the  above 
mathematically: 


ffi(y).n,  x  =  X\,  0  <y<L 

\F2(y).n,  x  =  x4,  0<y<L 


(38) 


where  Fi  (y)  and  F2(y)  represent  the  flow  field  solution  from  the  3D 
model.  Similar  to  the  3D  case,  the  outlet  is  specified  as  a  constant 
pressure  boundary: 


p  =  0{x2  <  x  <  x3 ,  y  =  0} 


(39) 


2.4.2.  Mass  transport  boundary  conditions 

Since  no  electrochemical  reactions  occur  outside  the  electrodes, 
the  concentration  is  assumed  constant  at  the  inlet  for  both  cathode 
and  anode  sides: 


r°  / . J  e  {IV,  V},  x  =  xi,  0  <  y  <  L 

1  j  \je  {//,///},  x  =  x4,  0  <y<L 


(40) 


The  high  purity  vanadium  solutions  in  4M  H2S04  used  in  [7] 
are  considered  for  which  c°  «  2  M,  %2M.  The  outlet  is  speci¬ 
fied  as  a  pure  convection  boundary  considering  the  relatively  large 
velocities  compared  to  other  regions  in  the  cell: 

dq 

^  =0{x2  <x<x3,  y  =  0}  (41) 

ay 


All  other  boundaries  are  considered  as  no-flux  or  insulation 
boundaries. 


2.4.3.  Current  conservation  boundary  conditions 

The  electric  current  enters  and  leaves  the  fuel  cell  via  the  two 
current  collectors  at  the  top  end  of  the  carbon  paper  electrodes.  For 
complete  cell  modeling  purposes,  the  anode  current  collector  is  set 
to  zero  potential  and  the  cathode  current  collector  is  maintained 
at  a  potential  which  is  varied  in  order  to  simulate  the  polarization 
curves  and  cell  performance: 


fVceih  *i<x<x2,  y  =  L 
|0,  x3  <  x  <  x4,  y  =  L 


(42) 


where  Vceu  is  the  cell  voltage.  Equivalently,  the  cathode  side  may  be 
set  to  a  current  density  boundary,  where  the  total  normal  current 
density  from  the  boundary  is  varied  to  simulate  cell  performance. 
For  half-cell  modeling,  the  electrolyte  potential  at  the  center  chan¬ 
nel  is  set  to  zero  and  the  potential  at  the  cathodic/anodic  current 
collector  is  swept: 


0s  =  Vceii{X\  <X  <x2,  y  =  L  or  x3<x<x4,  y  =  L]  (43) 
</>,  =  o{x=?^p,  0  <y  <  l)  (44) 

All  other  boundaries  are  set  to  insulation  boundaries. 


2.4.4.  Modeling  parameters 

Mass  transfer  and  electrochemical  kinetics  parameters  are  an 
important  part  of  the  numerical  model  in  the  current  work.  Identi¬ 
fication  of  relevant  and  accurate  values  for  these  is  closely  related 
to  model  accuracy.  A  thorough  search  of  the  literature  report¬ 
ing  experimental  measurements  was  conducted  to  ascertain  the 
parameter  values.  The  experimental  works  in  literature  [29,44-50] 
were  found  to  vary  based  on  conditions  such  as  vanadium  ion 
concentration,  acid  concentration,  electrode  material,  electrode 
surface  preparation  and  electrochemical  methods  used.  For  the 
vanadium  ions,  a  decreasing  trend  of  diffusion  coefficients  with 
increasing  concentration  of  the  vanadium  ion  and/or  acid  was 
found  [29,47,48].  The  viscosity,  however,  was  found  to  have  an 
increasing  trend  [51].  Variation  of  the  kinetic  rate  constants  with 
concentration,  electrode  material  and  electrode  surface  prepara¬ 
tion  methods  was  also  seen,  but  no  clear  trend  could  be  deduced 
[49,50]. 

Sum  et  al.  [44]  performed  cyclic  voltammetry  measurements  for 
the  V02+/VOj  couple  at  low  concentrations  of  vanadium  (0.055  M 
VOj  in  1 .8  M  H2S04).  Wen  et  al.  [47]  reported  diffusion  coefficients 
for  2  M  V02+  in  4  M  H2S04.  Direct  measurements  of  the  VOj  diffu¬ 
sion  coefficients  were  not  found  for  the  conditions  in  the  current 
work.  The  same  are  estimated  indirectly  using  the  corresponding 
Stokes  radius  [48].  The  relationship  between  Stokes  radius,  viscos¬ 
ity  and  diffusion  coefficient  can  be  expressed  as  [29]: 


D.  =  -W - 

1  6tt  fiRj 


(45) 


where  kb  is  the  Boltzmann  constant  and  Rj  is  the  Stokes  radius  of 
species  j.  For  the  V2+/V3+  couple,  Sum  et  al.  [45]  and  Yamamura 
et  al.  [49]  report  diffusion  coefficient  measurements  at  varied  solu¬ 
tion  conditions  and  electrode  materials.  However,  there  is  lack  of 
consistent  data  at  the  current  work’s  concentration  range.  There¬ 
fore,  the  diffusion  coefficients  of  V2+  and  V3+  are  calculated  using 
their  Stokes  radii.  The  viscosity  of  the  acidic  vanadium  solutions  is 
assumed  as  0.008  Pa  s  for  the  entire  fluid  domain  as  a  representative 
viscosity  for  all  four  vanadium  species.  The  charge  transfer  coeffi¬ 
cients  (a)  for  the  Butler-Volmer  equations  are  assumed  to  be  0.5 
for  both  half-cells.  A  summary  of  the  parameters  used  for  modeling 
is  provided  in  Tables  1-4. 
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Table  1 

Operating  conditions. 


Quantity 

Symbol 

Value 

Unit 

Source 

Operating  temperature 

T 

25 

°C 

- 

Volumetric  flow  rate 

V? 

in 

1-300 

1.667  x  10-11  to5x  10“9 

pi  min-1 

m3  s-1 

Measured  [7] 

Length 

L 

12  x  10“3 

m 

Measured 

Width 

W 

9  x  10-3 

m 

Measured 

Center  channel  width 

wc 

1  x  10“3 

m 

Measured 

Inlet  reservoir  width 

wr 

3x  10“3 

m 

Measured 

Inlet  channel  area 

O-in 

6x  10“7 

m2 

Measured 

Table  2 

Electrolyte  properties. 


Quantity 

Symbol 

Value 

Unit 

Source 

Density 

P 

1200 

kg  ITT3 

[29,51] 

Viscosity 

li 

0.008 

Pas 

[29,47,48] 

VOJ  diffusion  coefficient 

Dv 

0.364  xlO-10 

m2  s_1 

[29,47,48] 

V02+  diffusion  coefficient 

Div 

0.364  x  10-10 

m2  s_1 

[29,47,48] 

V3+  diffusion  coefficient 

Dm 

0.240  xlO”10 

m2  s_1 

Estimated 

V2+  diffusion  coefficient 

Du 

0.240  xlO”10 

m2  s_1 

Estimated 

Table  3 

Electrode  properties. 

Quantity 

Symbol 

Value 

Unit 

Source 

Active  length 

I 

KJ 

X 

o 

m 

Measured 

Width 

We 

1  X  10-3 

m 

Measured 

Thickness 

h 

300 x 10-6 

m 

[54] 

Porosity 

8 

0.78 

- 

[54] 

Fibre  diameter 

d 

7  x  10-6 

m 

[22] 

Pore  diameter 

dp 

28xl0-6 

m 

Estimated  [23] 

Permeability 

I< 

2.375  xlO”11 

m2 

[22] 

Electrode  platform  area 

A 

1.2  xl0“5 

m2 

Measured 

Specific  surface  area 

a 

8.333  xlO4 

m-1 

Estimated 

Electrical  conductivity 

of 

1.15  x  104 

Sm4 

[22] 

(in-plane) 

solution  procedure  highly  sensitive  to  the  initial  values  provided. 
Therefore,  in  order  to  achieve  stable  and  fast  convergence,  the  COM- 
SOL  segregated  solver  is  employed.  The  MUMPS  direct  linear  solver 
is  used  for  all  three  sub-models.  A  structured  mesh  with  boundary 
layers  is  used  for  both  the  3D  and  2D  models,  and  mesh  indepen¬ 
dence  and  convergence  are  ensured.  The  total  current  density  is 
obtained  by  numerical  integration  of  the  normal  electrode  current 
density  over  the  current  collector  boundary: 


(rs.h)dx 


(46) 


where  h  and  A  are  the  electrode  thickness  and  platform  area,  respec¬ 
tively  [7].  A  similar  numerical  integration  of  concentration  at  the 
outlet  boundary  is  performed  to  calculate  the  single  pass  fuel  uti¬ 
lization: 


(47) 


The  solution  procedure  is  repeated  over  a  selected  range  of 
applied  cell  or  electrode  potentials  to  obtain  complete  performance 
curves.  Typical  solution  times  for  the  3D  and  2D  models  are  around 
1 5  and  2  h,  respectively,  on  a  standard  PC. 


2.5.  Solution  procedure 

The  modeling  framework  described  in  the  previous  sections  is 
solved  based  on  conservation  of  mass,  momentum,  species  and 
charge.  Eqs.  (l)-(4)  are  first  solved  for  the  velocity  and  pressure 
distributions  in  the  microchannels  and  porous  media.  Eq.  (5)  is  then 
solved  for  the  concentration  of  each  of  the  four  vanadium  species 
together  with  the  charge  conservation  equation  (Eq.  (6))  to  find  the 
electrode  and  electrolyte  potentials  and  current  densities  at  a  given 
applied  cell  potential.  The  source  terms  for  species  and  charge  con¬ 
servation  are  coupled  through  Eqs.  (26)  and  (27).  This  coupling  also 
defines  the  relationship  between  surface  and  bulk  concentrations 
(Eqs.  (28)-(31)).  The  numerical  simulations  are  performed  using 
the  COMSOL  multiphysics  code  which  utilizes  a  Finite  Element 
Method  (FEM)  formulation.  A  quadratic  discretization  scheme  is 
used  for  all  the  mass  transport  and  electrochemical  variables.  The 
coupling  of  the  source  terms  between  the  mass  transport  and  elec¬ 
trochemistry  sub-models  introduces  non-linearity  and  makes  the 


3.  Results  and  discussion 

3.1.  Model  validation 

The  present  computational  model  is  developed  based  on  the 
design  of  the  microfluidic  fuel  cell  with  flow-through  porous  elec¬ 
trodes  previously  reported  by  our  group  [7].  To  characterize  the 
fundamental  physicochemical  phenomena  and  performance  of 
both  individual  half-cells  (anode  and  cathode)  as  well  as  the  com¬ 
bined  cell,  simulations  are  performed  for  the  entire  range  of  flow 
rates  considered  in  the  experiments,  i.e.  1-300  |xLmin_1.  The  first 
step  in  the  model  validation  process  is  to  validate  the  2D  sim¬ 
plification  introduced  in  Section  2.1  to  reduce  the  computational 
time  and  resources  required  for  the  simulations.  To  achieve  this, 
the  modeling  of  mass  transport  and  electrochemical  reactions  for 
a  2D  slice  are  compared  to  a  full  3D  model  of  the  cell.  Using  3D 
modeling  results  as  the  basis,  the  relative  deviation  between  the 
2D  and  3D  models,  averaged  over  the  modeling  domain,  is  found  to 
be  less  than  2%  when  compared  using  XY  distributions  of  electrode 


Table  4 

Electrochemical  parameters. 


Quantity 

Symbol 

Value 

Unit 

Source 

Standard  equilibrium  potential  (vs.  SCE) 

U? 

0.75 

V 

[30] 

Standard  equilibrium  potential  (vs.  SCE) 

u2° 

-0.496 

V 

[30] 

Standard  rate  constant 

fci 

6.8  xlO”7 

ms-1 

[49] 

Standard  rate  constant 

k2 

1.7  x  10“7 

ms4 

[45] 

Charge  transfer  coefficient  ( V02+  / VO  + -cathodic  reaction) 

au 

0.5 

- 

Assumed 

Charge  transfer  coefficient  ( V02+  / VO +  -anodic  reaction) 

Ofl, a 

0.5 

- 

Assumed 

Charge  transfer  coefficient  (V2+/V3+-cathodic  reaction) 

0(2, c 

0.5 

- 

Assumed 

Charge  transfer  coefficient  (V2+/V3+-anodic  reaction) 

0(2,  a 

0.5 

- 

Assumed 
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Fig.  2.  Comparison  of  simulated  polarization  curves  from  2D  and  3D  modeling 
methodology  at  a  moderate  flow  rate. 


current  density.  The  deviation  is  found  mainly  near  the  channel 
walls,  where  the  boundary  layer  predictions  of  the  3D  model  lead 
to  lower  current  densities  as  compared  to  the  2D  model.  The  two 
modeling  approaches  are  also  compared  using  concentration  and 
overpotential  distributions  which  show  the  same  relative  devia¬ 
tion  as  mentioned  above.  Polarization  curves  resulting  from  the 
two  approaches  are  presented  in  Fig.  2,  again  indicating  a  very 
close  agreement.  Other  flow  rate  cases  were  also  simulated  to 
confirm  these  trends,  thus  fully  justifying  the  2D  simplification 
of  the  proposed  model.  In  subsequent  sections,  results  based  on 
mass  transport  and  electrochemistry  are  discussed  using  a  2D  slice 
(Fig.  1(b))  of  the  electrodes  and  center  channel  at  the  symmetry 
z-plane,  whereas  the  entire  3D  geometry  is  used  to  discuss  fluid 
flow. 

Next,  we  seek  to  validate  the  numerical  model  with  experimen¬ 
tal  data  from  our  previous  study  [7].  For  a  relevant  comparison 
of  simulated  and  measured  data,  the  external  ohmic  resistance  of 
the  particular  fuel  cell  used  for  the  measurements  must  first  be 
assessed  and  added  to  the  modeling  parameters,  provided  it  is 
not  captured  within  the  modeling  domain  described  in  Fig.  1.  In 
the  measurements  reported  previously  [7],  the  electrical  contact 
between  the  carbon  electrodes  and  external  wires  was  made  using 
conductive  silver  epoxy,  which  is  expected  to  result  in  a  sizeable 
contact  resistance  owing  to  the  highly  porous  nature  of  the  elec¬ 
trode  surface.  A  simple  procedure  is  used  to  estimate  the  external 
resistance  from  the  total  ohmic  cell  resistance  measured  by  Elec¬ 
trochemical  Impedance  Spectroscopy  (EIS)  [7].  The  internal  ohmic 
resistance  is  first  calculated  based  on  electrode  and  electrolyte  con¬ 
ductivities  (Table  3  and  Eqs.  (33)  and  (34))  and  the  average  current 
path  inside  the  cell.  This  results  in  an  estimated  internal  ohmic 
resistance  (electrode  and  electrolyte)  of  5.9  £2.  The  total  ohmic 
resistance  of  the  cell,  including  external  contacts  and  wires,  was 
measured  to  be  31  £2  using  EIS  [7].  Using  this  value  with  the  internal 
resistance  estimated  above,  the  external  resistance  is  calculated  to 
be  25.1  £2.  Thus,  a  value  of  25.1  £2  is  added  externally  to  the  model  in 
order  to  compare  with  the  specific  experimental  setup  in  [7].  A  sim¬ 
ilar  procedure  can  be  adopted  to  estimate  the  external  resistance 
for  a  different  setup. 

After  adding  the  external  resistance,  the  simulated  results  can  be 
effectively  compared  to  the  measured  data.  It  is  important  to  note 
that  these  comparisons  are  based  on  zero  dimensional  data  such  as 
polarization  curves,  power  density  and  fuel  utilization  since  mea¬ 
sured  2D  distributions  of  concentration,  overpotential  and  current 
density  are  not  presently  available.  In  Fig.  3(a)  and  (b),  the  results  of 
half-cell  simulations  for  the  1 0  and  60  p,L  min-1  flow  rates  are  com¬ 
pared  to  the  corresponding  measured  polarization  data.  Complete 
cell  simulations  compared  to  measurements  for  both  low  (1  and 
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Fig.  3.  Comparison  of  half-cell  polarization  curves  predicted  by  the  model  with 
measurements  in  [7]:  (a)  10  p,Lmin_1  flow  rate  case;  and  (b)  60  p,Lmin_1  flow  rate 
case. 

10  pLmin-1 )  and  high  (60  and  300  pLmin-1 )  flow  rates  are  shown 
in  Fig.  4.  It  can  be  seen  that  good  agreement  is  achieved  for  both 
half-cell  and  fuel  cell  polarization  curves  across  the  whole  range  of 
flow  rates  considered.  The  agreement  is  seen  to  be  closer  at  higher 
flow  rates  as  the  sensitivity  of  the  model  to  mass  transport  parame¬ 
ters  is  lower  in  this  region.  The  average  relative  deviation  between 
model  and  measurements  across  all  flow  rates  considered  is  about 
1 0%.  The  maximum  deviation  is  found  to  be  1 9%  in  the  1 0  p,L  min-1 
case.  Figs.  5  and  6  show  comparisons  for  power  density  and  fuel 
utilization  with  measurements.  Good  agreement  is  observed,  with 
the  same  relative  deviation  as  mentioned  above.  The  model  gives  a 
good  estimate  of  the  peak  power  density  point  as  well  as  maximum 
fuel  utilization. 

Initial  measurements  based  on  new  cells  built  in  our  recently 
established  lab  at  Simon  Fraser  University  are  also  included  in  this 
study  for  comparison.  These  cells  are  comparable  in  design  to  the 
cell  built  at  the  University  of  Victoria  [7],  however,  the  channel  and 
electrode  thicknesses  are  now  reduced  to  180  pan  as  compared  to 
300  pm  in  the  previous  work.  As  shown  in  Fig.  7,  good  agreement 
is  again  observed  for  both  moderate  and  high  flow  rates.  The  model 
is  seen  to  follow  the  measured  trends  even  more  closely  for  this  set 
of  data.  The  relative  deviation  is  found  to  be  less  than  10%  for  the 
flow  rates  considered. 

It  is  important  to  gain  an  understanding  of  the  general  trends 
and  deviations  that  exist  between  model  and  measurements.  Some 
of  the  prominent  features  and  possible  reasons  for  the  deviations 
are: 

1.  Uncertainty  in  kinetic  and  mass  transport  parameters  obtained 
from  measured  values  reported  in  literature.  No  tuning  parame¬ 
ters  are  used  in  the  present  model  to  ensure  general  applicability 
with  respect  to  cell  chemistry  and  design. 
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Fig.  4.  Comparison  of  fuel  cell  polarization  curves  predicted  by  the  model  with 
measurements  in  [7]  at  four  different  flow  rates:  (a)  1,  10  (pLmin-1);  and  (b)  60, 
300  (pLmin-1). 
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Fig.  6.  Predicted  fuel  utilization  as  a  function  of  cell  voltage,  compared  to  mea¬ 
surements  in  [7]:  (a)  1,  10  (pLmin-1)  flow  rates;  and  (b)  60,  300  (pLmin-1)  flow 
rates. 
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Fig.  5.  Predicted  power  density  curves  compared  to  measurements  in  [7]:  (a)  1, 10 
(pi min-1)  flow  rates;  and  (b)  60,  300  (pLmin-1)  flow  rates. 


2.  Uncertainty  in  the  convective  mass  transfer  correlation  (Eq.  (21 )) 
of  the  order  of  30%  [38]. 

3.  The  model  assumes  a  uniform  pore  size  distribution  in  the  elec¬ 
trodes,  whereas  in  reality  the  carbon  paper  material  employed 
is  known  to  exhibit  a  high  degree  of  anisotropy  and  a  wide  range 
of  characteristic  pore  sizes  evident  from  microscope  images  and 
porosimetry  measurements  [39,52].  This  may  cause  channelling 
of  the  flow  through  predominantly  large  pores  and  reduced  per¬ 
formance  in  the  mass  transport  limited  regime. 

4.  The  small  deviation  at  low  current  densities  may  be  due  to  the 
assumption  of  single  electron  transfer  kinetics  for  the  vanadium 
reactions  [26]. 


Fig.  7.  Comparison  of  fuel  cell  polarization  curves  predicted  by  the  model  to  new 
measurements. 
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In  sum,  although  some  minor  deviations  are  found  between 
simulated  and  measured  data,  the  model  is  shown  to  predict  the 
trends  observed  in  measurements  within  a  reasonable  tolerance 
limit  and  provides  useful  estimates  of  half-cell  and  cell  polariza¬ 
tion,  power  density  and  fuel  utilization  over  the  entire  range  of 
current  densities  and  flow  rates  considered. 

While  the  fuel  cells  in  [7]  displayed  large  external  resistances, 
efforts  currently  underway  at  our  lab  indicate  that  the  exter¬ 
nal  resistance  can  be  brought  down  to  a  value  close  to  10  £2.  In 
subsequent  sections,  the  results  are  discussed  based  on  the  compu¬ 
tational  model  with  this  external  resistance  added.  This  is  to  ensure 
that  the  modeling  results  take  into  account  the  issue  of  large  contact 
resistances  which  are  present  in  these  cells  and  at  the  same  time 
maintain  generality  in  the  discussions  and  conclusions.  In  adapting 
these  general  results  to  a  specific  case,  it  is  advisable  to  measure 
external  resistances  for  the  specific  experimental  setup  and  adjust 
predicted  voltages  accordingly.  The  results  are  discussed  for  two 
main  operating  conditions:  (1)  maximum  efficiency  point;  and  (2) 
limiting  current  density  point  (mass  transport  limit). 

3.2.  Effect  of  flow  distribution 

The  3D  flow  distribution  in  the  electrodes  and  inlet/outlet 
channels  is  analyzed  for  uniformity.  Larger  cross-flow  velocity  mag¬ 
nitudes  (in  the  x  direction)  are  seen  in  the  portions  of  the  electrodes 
near  the  outlet.  These  magnitudes  are  found  to  be  up  to  2.5  times 
higher  than  those  near  the  current  collectors.  This  non-uniformity 
is  mainly  attributed  to  the  center  channel  which  is  narrower  in 
width  as  compared  to  the  inlet  channels  (Fig.  1(a))  thereby  con¬ 
tributing  a  larger  pressure  drop.  This  is  expected  to  impact  the 
concentration  and  potential  distributions  in  the  cell  due  to  differ¬ 
ences  in  rate  of  mass  transfer.  Flow  within  the  porous  electrodes 
is  effectively  uniform  in  the  z-direction  due  to  the  plug  flow  con¬ 
ditions  occurring  when  the  pore  size  is  significantly  smaller  than 
the  channel  height.  For  the  same  reason,  boundary  layers  within 
the  electrode  are  found  to  be  small  compared  to  those  in  the  free 
channel. 

The  parasitic  pumping  power  required  to  drive  the  electrolyte 
flow  can  be  analyzed  based  on  the  simulated  pressure  drops  in  the 
microfluidic  channels  and  porous  electrodes  of  the  cell  and  sub¬ 
tracted  from  the  total  power  produced.  The  total  pumping  power 
is  determined  by  the  product  of  the  pressure  drop  (from  inlet  to 
outlet)  and  the  volumetric  flow  rate  for  the  two  streams: 

Ppump  —  2A  pV?  (48) 

For  the  lOpuLmin-1  case,  the  pumping  power  is  found  to  be 
6.46  x  1 0-6  mW  and  for  the  300  fxL  min-1  case  it  is  5.82x1 0-3  mW. 
The  fraction  of  total  pumping  power  to  corresponding  power  out¬ 
put  of  the  cell  is  thus  less  than  0.1%  in  all  cases  considered,  and  can 
therefore  be  neglected. 

3.3.  Cell  performance  and  optimum  operating  point 

Fig.  8  presents  the  half-cell  and  complete  cell  polarization  curves 
obtained  by  the  model,  representing  the  theoretical  performance 
of  the  fuel  cell  for  the  case  of  10  £2  external  contact  or  wire  resis¬ 
tances.  It  can  be  seen  from  Fig.  8(a)  that  the  anode  is  the  limiting 
electrode  in  the  present  case.  The  reason  for  this  lies  in  the  lower 
diffusion  coefficient  values  of  the  anode  side  vanadium  species  as 
compared  to  the  cathode  which  means  that  transport  to  and  from 
the  electrode  is  slower  on  the  anode  side,  leading  to  a  lower  limit¬ 
ing  current  density  compared  to  the  cathode.  As  seen  in  Fig.  8(b), 
the  polarization  curves  for  all  flow  rates  considered  display  three 
distinct  regions:  an  activation  loss  region  at  low  current  densities;  a 
constant  ohmic  slope  at  mid-range  current  densities;  and  a  limiting 
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Fig.  8.  Predicted  theoretical  polarization  curves  with  external  resistance  of  10  £2\  (a) 
half-cell  polarization  curves  for  10  and  60  (p,Lmin_1 );  and  (b)  fuel  cell  polarization 
curves  for  1, 10,  60  and  300  (p,Lmin_1 )  flow  rates. 


current  drop  at  high  current  densities.  The  limiting  current  region  is 
generally  present  except  at  the  highest  flow  rate  (Fig.  8(b))  consid¬ 
ered  here  where  the  large  ohmic  polarization  causes  the  voltage  to 
drop  to  zero  before  the  limiting  current  density  is  attained.  The  lim¬ 
iting  current  value  follows  a  general  trend  with  flow  rate  consistent 
with  a  Levich  type  relation.  Due  to  the  high  in-plane  conductivity 
of  the  carbon  paper  electrodes  combined  with  a  high  concentra¬ 
tion  of  supporting  electrolyte,  the  cell  is  expected  to  have  a  low 
internal  resistance  as  reflected  in  the  gentle  slopes  in  the  ohmic 
region  of  the  curves.  These  findings  indicate  that  even  with  flow¬ 
through  porous  electrodes  the  mass  transport  to  the  active  sites  is 
still  the  main  constraint  to  cell  performance  (under  the  assumption 
that  the  external  contact  resistances  can  be  reduced  to  a  moderate 
value  (<10  £2)). 

Power  density  and  single  pass  fuel  cell  efficiency  are  two  of 
the  most  relevant  performance  metrics  for  microfluidic  fuel  cells 
with  flow-through  porous  electrodes,  and  for  most  applications  it 
is  desirable  to  achieve  concurrently  high  levels  of  these  two  char¬ 
acteristics.  In  the  present  flow-through  case,  the  power  density  is 
determined  by  the  product  of  current  density  and  cell  voltage,  while 
the  single  pass  fuel  cell  efficiency  is  determined  by  the  product  of 
current  and  voltage  efficiencies.  Provided  that  only  the  main  vana¬ 
dium  redox  reactions  (Eqs.  (7)  and  (8))  are  considered  in  the  model, 
the  current  efficiency  is  equal  to  the  fuel  utilization,  as  given  by: 


~  Sfuel  — 


1  cell A 


w, 


(49) 


The  voltage  efficiency  based  on  the  standard  cell  voltage  of 
U°  =  1.246  V  is  given  by: 


£v  = 


(VceA 

Wee,,) 


(50) 
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Fig.  9.  Concentration  distribution  of  V2+  (fuel)  across  the  center  channel  at  (a)  an  upstream  location  (y  =  L)  and  (b)  a  downstream  location  (y  =  0),  for  1, 10  and  300  (imLmirr1) 
flow  rates  at  OCV  conditions. 


The  single  pass  fuel  cell  efficiency  then  becomes: 


£ cell  —  x  — 


Ycell}cell^ 


(51) 


In  coherence  with  the  power  density,  the  single  pass  efficiency  is 
thus  proportional  to  the  product  of  current  density  and  cell  voltage. 
Consequently,  the  theoretical  maximum  power  density  and  single 
pass  efficiency  of  the  microfluidic  fuel  cell  always  occur  at  the  same 
operational  point  on  the  polarization  curve.  This  favourable  perfor¬ 
mance  characteristic  greatly  simplifies  cell  design  and  operation  as 
the  optimum  point  for  power  and  efficiency  is  shared.  In  the  present 
modeling  study,  the  maximum  cell  efficiency  is  found  to  be  77%  for 
1  puLmin-1  at  a  cell  voltage  of  1.1  V.  The  efficiency  is  found  to  peak 
at  high  cell  voltages  for  low  flow  rates  and  at  moderate  voltages 
for  high  flow  rates,  attributed  to  the  high  current  densities  and 
subsequently  high  ohmic  losses  at  high  flow  rates.  At  the  maxi¬ 
mum  efficiency  point,  the  theoretical  power  densities  increase  with 
flow  rate  from  26  rnW  cm  2  for  1  plmin-1  to  1 93  mW  cm  2  for 
300  pLmin-1.  The  high  power  densities  may  in  some  cases  justify 
operation  of  the  cell  at  high  flow  rates  despite  the  low  efficien¬ 
cies  (<10%),  e.g.,  for  applications  having  high  power  requirements 
for  short  durations,  while  a  notable  power  density  trade-off  is  evi¬ 
dent  for  operation  in  the  highly  efficient  low  flow  rate  regime.  In 
sum,  the  high  theoretical  levels  of  cell  efficiency  and  power  density 
predicted  by  the  model,  and  the  fact  that  these  two  characteristics 
share  the  same  optimum  operating  point  at  a  given  flow  rate,  is 
encouraging  towards  further  development  of  this  technology. 


3.4.  Secondary  reactions 

At  low  flow  rates  close  to  1  pLmin-1,  the  model  predicts  a 
lower  limiting  current  density  than  the  experimental  measure¬ 
ments  (Fig.  4(a))  which  tend  to  deviate  slightly  from  the  vertically 
drooping  polarization  curve  characteristic  of  the  mass  transfer 
limit.  In  addition,  the  measurements  at  1  plmin-1  show  fuel  uti¬ 
lization  greater  than  100%  (Fig.  6(a))  [7].  This  is  likely  attributed  to 
the  secondary  V3+/V 02+  redox  reactions  observed  in  measurements 
at  cell  voltages  below  0.4  V  [7].  The  modeling  results  are  used  to 
analyze  the  extent  and  effects  of  this  reaction.  Local  overpotentials 
at  the  cathode  and  anode  are  calculated  for  the  V3+/V 02+  couple 
(Eq.  (9))  as:  !]  =  &  -  fa-U,  where  U  =  0.07  V  vs.  SCE  [30].  At  low 
flow  rates  and  high  current  densities  (low  cell  voltages),  increased 
fuel  utilization  results  in  relatively  high  concentrations  of  product 
V02+  and  V3+  at  the  cathodic  and  anodic  side,  respectively.  Owing  to 
the  relatively  large  electrochemical  potential  gap  between  the  sec¬ 
ondary  reaction  and  the  two  primary  reactions  (0.57  and  0.68  V  for 
anode  and  cathode,  respectively),  sufficiently  high  overpotentials 


required  for  the  secondary  reaction  to  occur  are  only  found  for  cell 
voltages  below  ~0.6V.  At  this  cell  voltage,  the  concentration  of 
product  V02+  and  V3+  averaged  over  the  respective  electrode  area 
is  found  to  be  1.4  M,  which  is  high  enough  to  drive  the  reaction  for¬ 
ward  kinetically.  Thus,  both  thermodynamic  and  kinetic  conditions 
for  the  secondary  reaction  are  satisfied  when  the  cell  is  operated 
at  concurrently  low  levels  of  flow  rate  and  cell  voltage.  It  is  note¬ 
worthy,  however,  that  the  secondary  reaction  is  negligible  over  the 
majority  of  the  operational  range  of  the  cell,  including  its  optimum 
performance  point  at  peak  power  and  maximum  efficiency. 

3.5.  Concentration  distribution  and  crossover 


The  inter-diffusive  mixing  between  the  two  stratified  streams  in 
the  center  channel  is  analyzed  under  simulated  open  circuit  voltage 
(OCV)  conditions  at  the  lowest  flow  rate,  where  the  fuel  and  oxidant 
concentrations  are  the  highest  and  diffusive  mixing  is  most  severe. 
The  concentration  variation  of  V2+  (fuel)  across  the  width  of  the  cen¬ 
ter  channel  at  an  upstream  (y  =  L)  and  downstream  (y  =  0)  location 
is  plotted  in  Fig.  9(a)  and  (b),  respectively.  It  is  evident  that  even 
for  the  extreme  case  of  1  pLmin-1  flow  rate  at  OCV  conditions  the 
mixing  width  is  small  compared  to  the  width  of  the  center  channel. 
Further,  the  mixing  width  decreases  along  the  length  of  the  center 
channel  (Fig.  9(a)  and  (b)).  This  seemingly  counter-intuitive  but 
favourable  result  is  however  a  unique  feature  of  this  flow-through 
electrode  architecture.  In  contrast  to  the  planar  or  Y-channel  fuel 
cell  designs  [20,53],  the  fluid  velocity  increases  with  downstream 
position  along  the  center  channel,  i.e.,  towards  the  outlet,  due  to  the 
relatively  uniform  cross-flow  through  the  electrodes.  As  shown  in 
Fig.  9,  the  mixing  width  is  even  smaller  for  the  higher  flow  rates 
considered. 

On  computing  the  crossover  flux  from  one  electrode  to  the  other, 
a  second  favourable  feature  of  the  unique  flow-through  archi¬ 
tecture  is  discovered.  In  the  present  case,  the  flow  taking  place 
at  the  porous  electrode/center  channel  interface  is  orthogonal  to 
the  principal  axis  of  the  center  channel,  in  contrast  to  earlier  cell 
architectures  [9-12,14]  where  the  flow  is  predominantly  parallel 
to  the  principal  axis  in  flow-by  (or  flow-over)  mode.  Due  to  the 
flow-through  configuration,  any  crossover  flux  from  one  electrode 
towards  the  other  thus  needs  to  diffuse  upstream  into  the  oppo¬ 
site  electrolyte  flow  in  order  to  reach  the  opposite  electrode  and 
contribute  towards  secondary  reactions.  A  net  upstream  diffusive 
flux  can  only  occur  when  the  convective  flux  in  the  flow  direction 
is  lower  than  the  upstream  rate  of  diffusion.  The  resulting  Peclet 
number  for  crossover  flux  is  defined  based  on  channel  height: 


PpCo 
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Fig.  11.  Overpotential  (p)  distributions  along  the  electrode  length  for  cell  operation 
at  the  maximum  efficiency  point  (obtained  at  the  midpoint  of  the  electrode  width). 
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Fig.  10.  Concentration  profiles  of  V2+  (fuel)  in  the  symmetry  z-plane  at  selected  oper¬ 
ating  conditions  corresponding  to  the  maximum  efficiency  point  at  three  different 
flow  rates:  (a)  1  |jiLmin-1;  (b)  10  pi  min-1;  and  (c)  300  pLmin-1. 

where  v  ■  h  denotes  the  flow  velocity  normal  to  the  electrode 
boundary  (x  =  x3,  Fig.  1(a))  and  Dj  is  the  diffusion  coefficient  of  a 
given  vanadium  species.  Upstream  diffusion  is  expected  to  occur 
when  the  Peclet  number  defined  above  is  less  than  unity.  Simula¬ 
tions  of  the  fuel  cell  carried  out  at  very  low  flow  rates  (<1  pL  min-1 ) 
confirm  this  presumption.  The  diffusion  dominated  regime  is  found 
to  set  in  at  flow  rates  below  0.02  pL  min-1,  corresponding  to  a 
Peclet  number  of  0.1 .  At  flow  rates  above  this  value  the  net  flux  into 
the  electrodes  drops  sharply  with  increasing  flow  rate  due  to  the 
convection  dominated  mass  transport.  For  the  case  of  1  pLmin-1 
flow  rate  at  OCV  conditions,  the  net  flux  of  primary  reactant  into 
the  electrode  is  computed  to  be  three  orders  of  magnitude  greater 
than  the  crossover  flux.  Consequently,  crossover  flux  and  associ¬ 
ated  parasitic  secondary  current  can  be  safely  neglected  for  the  cell 
geometry  and  flow  rates  considered  in  this  study. 

The  concentration  plot  for  1  pLmin-1  flow  rate  obtained  at  the 
maximum  efficiency  point  (Fig.  1 0(a))  shows  that  most  of  the  fuel  is 
consumed  before  the  cross-flow  exits  the  electrode  into  the  center 
channel.  Furthermore,  a  major  contribution  for  the  reaction  appears 
to  occur  in  the  entrance  portion  of  the  electrode  where  the  reactant 


concentration  is  high.  The  concentration  profile  for  this  case  also 
reveals  that  the  electrode  is  underutilized  for  this  flow  rate  since 
half  of  the  electrode  width  contributes  to  90%  of  the  conversion. 
In  Fig.  10(b),  decreased  fuel  conversion  at  increased  flow  rates  is 
observed  for  the  1 0  pL  min-1  case.  As  flow  rates  increase,  the  higher 
current  densities  are  expected  to  make  ohmic  polarization  a  more 
dominant  factor.  In  Fig.  10(c)  for  the  300  pLmin-1  flow  rate,  the 
fuel  conversion  is  greater  near  the  current  collector  as  compared 
to  the  portion  of  the  electrode  near  the  outlet,  owing  to  the  ohmic 
potential  drop  along  the  electrode  length.  In  addition,  the  overall 
consumption  of  fuel  is  notably  smaller  at  this  flow  rate. 

3.6.  Electrode  performance  at  maximum  efficiency 

The  simulated  electrode  performance  characteristics  are  ana¬ 
lyzed  in  more  detail  for  operation  at  the  maximum  efficiency  point. 
The  local  overpotential  is  plotted  along  the  electrode  length  for 
a  low  (1  pLmin-1)  and  high  (300  pLmin-1)  flow  rate  in  Fig.  11. 
Increase  in  ohmic  polarization  effects  with  increased  flow  rate  is 
evident  due  to  larger  current  densities  generated  at  higher  flow 
rates.  At  low  and  moderate  flow  rates  the  distributions  have  a 
negligible  gradient  along  the  electrode  length,  whereas  they  show 
clear  variation  along  the  length  for  high  flow  rates  (Fig.  11).  For  the 
300  pL  min-1  case,  the  difference  in  overpotential  between  the  cur¬ 
rent  collector  and  the  opposite  end  of  the  electrode  is  on  the  order 
of  10-100  mV,  which  causes  the  trend  in  fuel  conversion  explained 
above  (Section  3.5).  The  local  overpotential  in  the  electrode  is  one  of 
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Fig.  12.  Normalized  current  density  distributions  in  the  symmetry  z-plane  at  the  maximum  efficiency  point:  (a)  1  p.Lmin-1;  (b)  lOjjiLmin-1;  and  (c)  300  jjiLmin-1. 
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Fig.  13.  Normalized  current  density  distributions  in  the  symmetry  z-plane  at  the  limiting  current  density:  (a)  1  (jiLmin-1 ;  (b)  10  jjiLmirr1 ;  and  (c)  300  jjiLmin-1. 


the  parameters  governing  the  reaction  rate  or  local  current  density 
at  that  point,  which  is  an  indicator  of  electrode  performance.  Other 
factors  are  the  local  velocity  magnitude,  concentration  difference 
between  and  bulk  and  surface,  as  well  as  the  reaction  rate  con¬ 
stants.  To  understand  the  combined  effect  of  these,  the  normalized 
current  density,  or  source  term  (Si  and  S2),  is  plotted  for  the  three 
different  flow  rates  in  Fig.  12.  This  current  density  is  normalized 
by  dividing  by  its  average  value  over  the  entire  electrode  platform 
area.  For  the  cathode  and  anode  side  it  is  given  as: 


jn=  lflhl 

1  Jf(\aU\/Lwe)dA 

,n  !f”2l 

2  fJ(\ai2\/Lwe)-dA 


(53) 

(54) 


As  per  the  normalization  method,  local  values  greater  than  the  aver¬ 
age  translate  to  normalized  values  greater  than  unity  and  vice  versa. 
For  ideal,  uniform  electrode  performance,  the  normalized  values 
should  be  close  to  unity  everywhere  in  the  electrode.  In  Fig.  12(a) 
and  (b),  representing  low  and  moderate  flow  rates,  a  visible  band 
of  high  local  current  density  values  is  observed  at  the  entrance  of 
the  electrode.  The  contours  at  low  flow  rates  (hence  low  current 
densities)  also  indicate  that  the  concentration  difference  between 
bulk  and  surface  is  the  main  factor  that  influences  the  variations  in 
local  current  density  along  the  electrode.  At  high  flow  rates,  how¬ 
ever,  the  overpotential  is  found  to  play  a  more  important  role,  with 
the  reaction  current  density  having  a  larger  value  near  the  current 
collectors  (Fig.  12(c)). 


3. 7.  Electrode  performance  at  limiting  current 

It  is  interesting  to  see  if  the  trends  discussed  above  at  maximum 
efficiency  also  hold  when  one  or  both  of  the  electrodes  is  produc¬ 
ing  the  maximum  or  limiting  current  density  for  a  given  flow  rate. 
Once  again  the  normalized  current  density  is  utilized  as  a  cumula¬ 
tive  indicator  of  the  various  factors  mentioned  above.  The  results 
in  Fig.  13(a)  and  (b)  are  similar  to  the  maximum  efficiency  plots  at 
low  and  moderate  flow  rates,  however,  a  reversal  in  trends  is  found 
at  high  flow  rates.  At  low  and  moderate  flow  rates,  the  performance 
is  predominantly  mass  transport  controlled  at  both  maximum  effi¬ 
ciency  and  limiting  current  density  points,  and  exhibits  a  slightly 
greater  current  contribution  from  the  region  of  the  electrodes  with 
greater  local  velocities.  At  high  flow  rates,  it  was  discussed  above 


that  the  overpotential  effects  dominate  hence  leading  to  greater 
contribution  from  the  regions  near  the  current  collector.  The  oppo¬ 
site  effect  is  found  at  the  limiting  current  density  point,  where  the 
regions  experiencing  greater  flow  velocities  (Section  3.2)  display 
larger  reaction  source  terms  despite  the  relatively  low  overpoten¬ 
tial  magnitudes.  This  reversed  trend  is  only  observed  in  the  anode, 
which  is  presently  the  limiting  electrode  of  the  cell.  It  is  deduced 
from  this  observation  that  local  velocities  play  a  significant  role 
in  determining  the  electrode  performance  at  the  limiting  current 
density  point. 

4.  Conclusions 

A  detailed  computational  model  of  a  microfluidic  fuel  cell 
with  flow-through  porous  electrodes  is  developed  and  validated. 
The  present  model  is  the  first  comprehensive,  cell-level  model 
of  this  fuel  cell  architecture.  The  three  major  coupled  phenom¬ 
ena  occurring  in  the  cell,  namely  the  fluid  flow,  mass  transfer  and 
electrochemical  kinetics,  are  solved  in  3D  using  a  computational 
multiphysics  solver.  Performance  measurements  carried  out  on 
these  cells  are  used  to  validate  the  resulting  model.  The  model¬ 
ing  results  are  found  to  be  in  good  agreement  with  experimental 
data  over  three  orders  of  magnitude  of  flow  rates  and  current  den¬ 
sities  spanning  the  complete  operational  range  of  these  devices. 
The  results  are  employed  to  verify  and  quantify  the  relatively  high 
external  contact  resistance  of  the  present  prototype  cells.  With 
a  representative  external  resistance  of  10  £2,  the  model  predicts 
that  theoretical  efficiencies  up  to  77%  and  power  densities  up  to 
193  mW  cm-2  are  feasible  at  low  and  high  flow  rates,  respectively. 
Interestingly,  the  maximum  efficiency  and  peak  power  of  this  tech¬ 
nology  are  found  to  share  the  same  optimum  operational  point  on 
the  polarization  curve  at  a  given  flow  rate.  Furthermore,  the  model¬ 
ing  results  indicate  that  secondary  reactions  and  reactant  crossover 
can  be  neglected  for  the  cell  architecture  and  operating  conditions 
considered  in  this  study. 

An  initial  investigation  is  carried  out  using  the  validated  model 
to  identify  variations  in  performance  within  the  porous  electrodes. 
The  analysis  based  on  simulated  spatial  distributions  of  concen¬ 
tration,  overpotential  and  current  density  reveals  some  interesting 
and  useful  trends.  Local  concentration  variations  determine  elec¬ 
trode  performance  for  low  and  moderate  flow  rates,  indicating 
mass  transfer  controlled  operation.  At  high  flow  rates,  both  mass 
transfer  and  ohmic  effects  are  significant,  and  highly  coupled  trends 


I<.  Deepak  et  al.  /  Journal  of  Power  Sources  196(2011)  10019-10031 


10031 


are  observed.  Cell  design  attributes  such  as  optimal  width  or  active 
area  of  the  electrode  to  achieve  high  fuel  utilization  (>90%)  are 
thus  highly  dependent  on  flow  rate  even  for  porous  flow-through 
architectures. 

Based  on  our  preliminary  investigation,  we  can  conclude  that 
performance  increases  are  possible  by  using  modeling  results  to 
drive  the  fuel  cell  design  cycle.  Some  possible  areas  of  improvement 
include  modifying  the  electrode  and  channel  geometries  to  reduce 
ohmic  polarization  effects  and  optimizing  the  electrolyte  flow  path 
based  on  flow  rate  and  required  power  density.  These  increases  are 
over  and  above  those  achieved  by  solving  the  contact  resistance 
problem  which  is  presently  the  main  bottleneck  for  performance 
improvement.  The  model  is  thus  useful  for  both  performance 
predictions  and  scientific  studies  to  gain  a  clear  understanding 
between  the  interplay  of  various  physicochemical  mechanisms  and 
their  effect  on  cell  performance.  The  model  is  also  anticipated  to 
become  a  useful  tool  for  design  and  analysis  of  other  fuel  cells 
and  electrochemical  sensors  which  incorporate  microfluidic  flow¬ 
through  porous  electrodes.  A  detailed  parametric  analysis  of  the 
flow-through  microfluidic  fuel  cell  based  on  the  model  developed 
here  will  be  presented  in  a  forthcoming  work. 
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